
conf = omni_load_conf()
s = omni_read_cat(conf.cat,ncat)

restore,'./local/'+conf.survey+'_velocities.sav',/ver
IF ~exist(pvec) THEN restore,'./local/'+conf.survey+'_pvec.sav',/ver

;; GRS overlap stats
ind = where(v.l GE 18 AND v.l LE 56, nv)
print,'Objects in the GRS overlap: ',nv
v = v[ind]
ind = where(v.vlsr GE -200, nvv)
print,'Valid VLSR in GRS overlap: ',nvv
v = v[ind]
ind = where(v.vlsr LE 0, nog)
print,'Objects in GRS overlap beyond the Solar circle: ',nog

;; GRS landlocked stats
restore,'./local/'+conf.survey+'_velocities.sav'
ind = where(v.grs.flag GT 3, nll)
print,'N GRS that are landlocked: ',nll

;; Velocity Detection stats
ind = where(v.vlsr GE -200, nval)
print,'N valid VLSR from dense gas: ',nval
ind = where((v.mol[0].multiv EQ 0 AND v.mol[0].nspec GE 1) OR $
            (v.mol[1].multiv EQ 0 AND v.mol[1].nspec GE 1), nhht)
print,'N sources with HHT detection: ',nhht
ind = where(v.mol[3].multiv EQ 0 AND v.mol[3].nspec GE 1, nnh3)
print,'N sources with NH3 detection: ',nnh3

;; GRS detection stats
f = v.grs.flag
ygrs = (f EQ 1) OR (f EQ 2) OR (f EQ 5) OR (f EQ 6)
ind = where(ygrs, ngrs)
print,'N sources with valid GRS MES: ',ngrs
ind = where(constrain.grsmatch, ngrsm)
print,'N sources with GRSMATCH: ',ngrsm
ind = where(v.vlsr GE -200 OR ygrs, nspec)
print,'N sources with EITHER dense gas OR 13CO: ',nspec
print,'Fraction with spectroscopic info: ',float(nspec)/ncat

;; DPDF Numbers:
print,'N sources with KDIST:    ',fix(total(constrain.kdist)),$
      fix(total(constrain.kdist AND constrain.post))
print,'N sources with GRSMATCH: ',fix(total(constrain.grsmatch)),$
      fix(total(constrain.grsmatch AND constrain.post))
print,'N sources with H2:       ',fix(total(constrain.h2)),$
      fix(total(constrain.h2 AND constrain.post))
print,'N sources with EMAF:     ',fix(total(constrain.emaf)),$
      fix(total(constrain.emaf AND constrain.post))
print,'N sources with PARALLAX: ',fix(total(constrain.parallax)),$
      fix(total(constrain.parallax AND constrain.post))
print,'N sources with HRDS:     ',fix(total(constrain.hrds)),$
      fix(total(constrain.hrds AND constrain.post))
print,'N sources with KNOWND:   ',fix(total(constrain.knownd)),$
      fix(total(constrain.knownd AND constrain.post))
print,'N sources with POST:     ',fix(total(constrain.post)),$
      fix(total(constrain.post AND constrain.post))

nwc  = fix(total(constrain.post))
ksam = constrain.kdist OR constrain.grsmatch OR $
                constrain.parallax OR constrain.knownd
nkin = fix(total(ksam))
print,'N sources with kinematics: ',nkin
print,'Fraction of kinematics with w.c.: ',float(nwc)/nkin

;; KDARs
near = fix(total(pvec.stat.kdar EQ 'N'))
far  = fix(total(pvec.stat.kdar EQ 'F'))
tan  = fix(total(pvec.stat.kdar EQ 'T'))
out  = fix(total(pvec.stat.kdar EQ 'O'))
excl = fix(total(pvec.stat.kdar EQ 'X'))
unc  = fix(total(pvec.stat.kdar EQ 'U'))
print,'N near: ',near,float(near)/nwc
print,'N far:  ',far, float(far)/nwc
print,'N tan:  ',tan, float(tan)/nwc
print,'N out:  ',out ,float(out)/nwc
print,'N excl: ',excl,float(excl)/nwc
print,'N unc : ',unc ,float(unc)/nwc

ind = where(constrain.post AND pvec.glon LE 90., nig)
print,'N w.c. in the inner Galaxy: ',nig

ind = where(constrain.post AND pvec.stat.duse[0]/1.d3 LE 5.5, n55)
print,'N w.c. nearer than 5.5 kpc: ',n55,float(n55)/nwc
ind = where(constrain.post AND pvec.stat.duse[0]/1.d3 LE 6.0, n60)
print,'N w.c. nearer than 6.0 kpc: ',n60,float(n60)/nwc


hl = abs(constrain.glat) GE 0.4
ihl = where( hl, nhl)

;; H2 Prior
nh2 = fix(total(constrain.h2))
ind = where(constrain.h2 AND (constrain.kdist OR constrain.grsmatch OR $
                              constrain.knownd OR constrain.parallax), nh2k)
print,'H2 Prior (total / kin): ',nh2,nh2k,float(nh2k)/nh2
ind = where(constrain.h2 AND (constrain.kdist OR constrain.grsmatch OR $
                              constrain.knownd OR constrain.parallax) AND $
            constrain.post, nh2wc)
print,'H2 Prior w.c.: ',nh2wc,float(nh2wc)/nh2k,float(nh2wc)/nh2
ind = where(hl AND constrain.h2,nh2hl)
ind = where(hl AND constrain.h2 AND constrain.post, nh2hlwc)
print,'High-Latitude: ',nhl,nh2hl,nh2hlwc,double(nh2hlwc)/nh2hl


END
